home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / clang / jcool01.zip / COMPLEX.C < prev    next >
C/C++ Source or Header  |  1992-09-28  |  15KB  |  435 lines
  1. //
  2. // Copyright (C) 1991 Texas Instruments Incorporated.
  3. // Copyright (C) 1992 General Electric Company.
  4. //
  5. // Permission is granted to any individual or institution to use, copy, modify,
  6. // and distribute this software, provided that this complete copyright and
  7. // permission notice is maintained, intact, in all copies and supporting
  8. // documentation.
  9. //
  10. // Texas Instruments Incorporated and General Electric Company
  11. // provides this software "as is" without express or implied warranty.
  12. //
  13. //
  14. // Created: MBN 11/01/89 -- Initial implementation
  15. // Updated: MBN 03/04/90 -- Added exception for DIVIDE_BY_ZERO
  16. // Updated: MJF 03/12/90 -- Added group names to RAISE
  17. // Updated: MJF 07/31/90 -- Added terse print
  18. // Updated: DLS 03/22/91 -- New lite version
  19. // Updated: VDN 06/29/92 -- roots of real polynomial, degree <= 4.
  20. //
  21. // The Complex  class implements  Complex  numbers  and arithmetic.   A Complex
  22. // object has the same  precision and range of  values  as  the system built-in
  23. // type double.  Implicit conversion  to  the  system defined types short, int,
  24. // long, float,    and  double   is supported by   overloaded  operator  member
  25. // functions.  Although   the  Complex class  makes   judicous use   of  inline
  26. // functions and deals only with floating point values, the user is warned that
  27. // the Complex double arithmetic class is still  slower than  the built-in real
  28. // data types.
  29. //
  30. // Find roots of a real polynomial in a single variable, with degree <=4.
  31. // Reference: Winston, P.H, Horn, B.K.P. (1984) "Lisp", Addison-Wesley.
  32.  
  33. #ifndef COMPLEXH                // If no Complex class
  34. #include <cool/Complex.h>            // Include class definition
  35. #endif
  36.  
  37. // minus_infinity -- Raise Error exception for negative infinity
  38. // Input:            Character string of derived class and function
  39. // Output:           None
  40.  
  41. void CoolComplex::minus_infinity (const char* name) const {
  42.   //RAISE (Error, SYM(CoolComplex), SYM(Minus_Infinity),
  43.   printf ("CoolComplex::%s: Operation results in negative infinity value",
  44.       name);
  45.   abort ();
  46. }
  47.  
  48.  
  49. // plus_infinity -- Raise Error exception for positive infinity
  50. // Input:           Character string of derived class and function
  51. // Output:          None
  52.  
  53. void CoolComplex::plus_infinity (const char* name) const {
  54.   //RAISE (Error, SYM(CoolComplex), SYM(Plus_Infinity),
  55.   printf ("CoolComplex::%s: Operation results in positive infinity value",
  56.       name);
  57.   abort ();
  58. }
  59.  
  60.  
  61. // overflow -- Raise Error exception for overflow occuring during conversion
  62. // Input:      Character string of derived class and function
  63. // Output:     None
  64.  
  65. void CoolComplex::overflow (const char* name) const {
  66.   //RAISE (Error, SYM(CoolComplex), SYM(Overflow),
  67.   printf ("CoolComplex::%s: Overflow occured during type conversion", name);
  68.   abort ();
  69. }
  70.  
  71.  
  72. // underflow -- Raise Error exception for underflow occuring during conversion
  73. // Input:       Character string of derived class name and function
  74. // Output:      None
  75.  
  76. void CoolComplex::underflow (const char* name) const {
  77.   //RAISE (Error, SYM(CoolComplex), SYM(Underflow),
  78.   printf ("CoolComplex::%s: Underflow occured during type conversion",name);
  79.   abort ();
  80. }
  81.  
  82.  
  83. // no_conversion -- Raise Error exception for no conversion from CoolComplex
  84. // Input:           Character string of derived class name and function
  85. // Output:          None
  86.  
  87. void CoolComplex::no_conversion (const char* name) const {
  88.   //RAISE (Error, SYM(CoolComplex), SYM(No_Conversion),
  89.   printf ("CoolComplex::%s: No conversion from CoolComplex", name);
  90.   abort ();
  91. }
  92.  
  93.  
  94. // divide_by_zero -- Raise Error exception for divide by zero
  95. // Input:            Character string of derived class and function
  96. // Output:           None
  97.  
  98. void CoolComplex::divide_by_zero (const char* name) const {
  99.   //RAISE (Error, SYM(CoolComplex), SYM(Divide_By_Zero),
  100.   printf ("CoolComplex::%s: Divide by zero", name);
  101.   abort ();
  102. }
  103.  
  104.  
  105. // operator= -- Overload the assignment operator for the CoolComplex class
  106. // Input:       Reference to CoolComplex object
  107. // Output:      Reference to updated CoolComplex object
  108.  
  109. CoolComplex& CoolComplex::operator= (const CoolComplex& c) {
  110.   this->r = c.r;                // Copy real part
  111.   this->i = c.i;                // Copy imaginary part
  112.   this->state = c.state;            // Copy state
  113.   return *this;                    // Return reference
  114. }
  115.  
  116.  
  117. // operator short -- Provide implicit conversion from CoolComplex to short
  118. // Input:            None
  119. // Output:           Converted number
  120.  
  121. CoolComplex::operator short () {
  122.   if (this->i != 0.0) {                // If there is an i-part
  123.     this->state = N_NO_CONVERSION;        // Indicate exception case
  124.     this->no_conversion ("operator short");    // And raise exception
  125.   }
  126.   else if (this->r > MAXSHORT) {        // If there is an overflow
  127.     this->state = N_OVERFLOW;            // Indicate exception case
  128.     this->overflow ("operator short");        // And raise exception
  129.   }
  130.   else if (this->r < MINSHORT) {        // If there is an underflow
  131.     this->state = N_UNDERFLOW;            // Indicate exception case
  132.     this->underflow ("operator short");        // And raise exception
  133.   }
  134.   return short(this->r);            // Return converted number
  135. }
  136.  
  137.  
  138. // operator int -- Provide implicit conversion from CoolComplex to int
  139. // Input:          None
  140. // Output:         Converted number
  141.  
  142. CoolComplex::operator int () {
  143.   if (this->i != 0.0) {                // If there is an i-part
  144.     this->state = N_NO_CONVERSION;        // Indicate exception case
  145.     this->no_conversion ("operator int");    // And raise exception
  146.   }
  147.   else if (this->r > MAXINT) {            // If there is an overflow
  148.     this->state = N_OVERFLOW;            // Indicate exception case
  149.     this->overflow ("operator int");        // And raise exception
  150.   }
  151.   else if (this->r < MININT) {            // If there is an underflow
  152.     this->state = N_UNDERFLOW;            // Indicate exception case
  153.     this->underflow ("operator int");        // And raise exception
  154.   }
  155.   return int (this->r);                // Return converted number
  156. }
  157.  
  158.  
  159. // operator long -- Provide implicit conversion from CoolComplex to long
  160. // Input:           None
  161. // Output:          Converted number
  162.  
  163. CoolComplex::operator long () {
  164.   if (this->i != 0.0) {                // If there is an i-part
  165.     this->state = N_NO_CONVERSION;        // Indicate exception case
  166.     this->no_conversion ("operator long");    // And raise exception
  167.   }
  168.   else if (this->r > MAXLONG) {            // If there is an overflow
  169.     this->state = N_OVERFLOW;            // Indicate exception case
  170.     this->overflow ("operator long");        // And raise exception
  171.   }
  172.   else if (this->r < MINLONG) {            // If there is an underflow
  173.     this->state = N_UNDERFLOW;            // Indicate exception case
  174.     this->underflow ("operator long");        // And raise exception
  175.   }
  176.   return long (this->r);            // Return converted number
  177. }
  178.  
  179.  
  180. // operator float -- Provide implicit conversion from CoolComplex to float
  181. // Input:            None
  182. // Output:           Converted number
  183.  
  184. CoolComplex::operator float () {
  185.   if (this->i != 0.0) {                // If there is an i-part
  186.     this->state = N_NO_CONVERSION;        // Indicate exception case
  187.     this->no_conversion ("operator float");    // And raise exception
  188.   }
  189.   else if (this->r > MAXFLOAT) {        // If there is an overflow
  190.     this->state = N_OVERFLOW;            // Indicate exception case
  191.     this->overflow ("operator float");        // And raise exception
  192.   }
  193.   else if (this->r < 0.0 && (-this->r) < MINFLOAT) { // Is there an underflow?
  194.     this->state = N_UNDERFLOW;            // Indicate exception case
  195.     this->underflow ("operator float");        // And raise exception
  196.   }
  197.   return float (this->r);            // Return converted number
  198. }
  199.  
  200.  
  201. // operator double -- Provide implicit conversion from CoolComplex to double
  202. // Input:             None
  203. // Output:            Converted number
  204.  
  205. CoolComplex::operator double () {
  206.   if (this->i != 0.0) {                // If there is an i-part
  207.     this->state = N_NO_CONVERSION;        // Indicate exception case
  208.     this->no_conversion ("operator double");    // And raise exception
  209.   }
  210.   return double (this->r);            // Return converted number
  211. }
  212.  
  213. // operator/= -- Overload the division/assign operator for the CoolComplex class
  214. // Input:       Reference to complex object
  215. // Output:      Reference to new complex object
  216.  
  217. CoolComplex& CoolComplex::operator/= (const CoolComplex& c2) {
  218.   if (c2.r == 0.0 && c2.i == 0.0)        // If both num and den zero
  219.     this->divide_by_zero ("operator/");        // Raise exception
  220.   if (c2.r == 0.0)                // If zero real part
  221.     if (this->r < 0.0 && this->i >= 0.0)    // If negative complex
  222.       this->minus_infinity ("operator/");    // Raise exception
  223.     else
  224.       this->plus_infinity ("operator/");    // Raise exception
  225.   if (c2.i == 0.0)                // If zero real part
  226.     if (this->i < 0.0 && this->r >= 0.0)    // If negative complex
  227.       this->minus_infinity ("operator/");    // Raise exception
  228.     else
  229.       this->plus_infinity ("operator/");    // Raise exception
  230.   double normalize = (c2.r * c2.r) + (c2.i * c2.i);
  231.   double new_r = (this->r * c2.r) + (this->i * c2.i); // multiply with conjugate
  232.   double new_i = (this->i * c2.r) - (this->r * c2.i);
  233.   this->r = new_r / normalize;
  234.   this->i = new_i / normalize;
  235.   return *this;
  236. }
  237.  
  238.  
  239. // operator<< -- Overload the output operator for a reference to a CoolComplex
  240. // Input:        Ostream reference, reference to a CoolComplex object
  241. // Output:       Ostream reference
  242.  
  243. ostream& operator<< (ostream& os, const CoolComplex& c) {
  244.   os << "(" << c.r << "," << c.i << ")";
  245.   return os;
  246. }
  247.  
  248.  
  249. // print --  terse print function for CoolComplex
  250. // Input:    reference to output stream
  251. // Output:   none
  252.  
  253. void CoolComplex::print(ostream& os) {
  254.   os << "                    /* CoolComplex %lx */" << (unsigned long) this;
  255. }
  256.  
  257.  
  258. // curt -- Cubic root of a double
  259.  
  260. double curt (double d) {
  261.   if (d == 0.0) 
  262.     return 0.0;
  263.   else if (d < 0) 
  264.     return -pow(-d, 1.0/3.0);
  265.   else
  266.     return pow(d, 1.0/3.0);
  267. }
  268.  
  269. // roots of linear polynomial: a*x + b = 0.
  270.  
  271. int CoolComplex::roots_of_linear (const double& a, const double& b, 
  272.                   CoolComplex& r) {
  273.   if (a == 0) {
  274.     if (b == 0) {
  275.       cout << "Homogenous has infinite number of roots." << endl;
  276.       r = 0;                    // least norm solution
  277.       return 1;
  278.     } else {
  279.       cout << "Inconsistent equation. No root." << endl;
  280.       return 0;
  281.     }
  282.   } else {
  283.     r = -b / a;                    // normal case, degree=1
  284.     return 1;                    
  285.   }
  286. }
  287.  
  288. // roots of quadratic -- a*x^2 + b*x + c = 0.
  289. //         Return two roots, largest magnitude first.
  290.  
  291. int CoolComplex::roots_of_quadratic (const double& a, const double& b, 
  292.                      const double& c, 
  293.                      CoolComplex& r1, CoolComplex& r2) {
  294.   if (a < 0) {                    // a is make positive, try again.
  295.     return roots_of_quadratic(-a, -b, -c, r1, r2);
  296.   } else if (a == 0) {                // special cases
  297.     return roots_of_linear(b, c, r1);
  298.   } else if (c == 0) {
  299.     r2 = 0;
  300.     return roots_of_linear(a, b, r1) + 1;
  301.   } else {                    // normal case
  302.     double discriminant = (b * b) - (4.0 * a * c);
  303.     double a2 = 2 * a;
  304.     if (discriminant < 0) {
  305.       double real = -b / a2;
  306.       double imag = sqrt(-discriminant) / a2;
  307.       r1 = CoolComplex(real, imag);        // two complex roots
  308.       r2 = CoolComplex(real, -imag);
  309.     } else if (discriminant == 0) {        // one double root
  310.       r1 = r2 = -b / a2;
  311.     } else {                    // two real roots
  312.       double n;
  313.       if (b < 0) n = sqrt(discriminant) - b;
  314.       else       n = -(sqrt(discriminant) + b);
  315.       r1 = n / a2;                // root with largest 
  316.       r2 = (2 * c) / n;                // magnitude first.
  317.     }
  318.     return 2;                    // number of roots=2
  319.   }
  320. }
  321.  
  322. // roots of cubic -- a*x^3 + b*x^2 + c*x + d = 0.
  323. //         Return three roots, largest magnitude first.
  324.  
  325.  
  326. int CoolComplex::roots_of_cubic (const double& a, const double& b, 
  327.                  const double& c, const double& d,
  328.                  CoolComplex& r1, CoolComplex& r2, 
  329.                  CoolComplex& r3) {
  330.   if (a < 0) {
  331.     return roots_of_cubic(-a, -b, -c, -d, r1, r2, r3);
  332.   } else if (a == 0) {                // special cases
  333.     return roots_of_quadratic(b, c, d, r1, r2);
  334.   } else if (d == 0) {
  335.     r3 = 0;
  336.     return roots_of_quadratic(a, b, c, r1, r2) + 1;
  337.   } else {                    // normal case
  338.     CoolComplex rs1, rs2;            // roots of resolvent
  339.     roots_of_quadratic(1, 
  340.                ((2*b*b*b) + (9 * a * ((3*a*d) - (b*c)))),
  341.                pow((b*b) - (3*a*c), 3),
  342.                rs1, rs2);
  343.     if (rs1.imaginary() == 0) {            // resolvent roots are real
  344.       double r = curt(rs1.real());        // find cube roots of resolvents
  345.       double s = curt(rs2.real());
  346.       double a3 = 3 * a;
  347.       r1 = (r + s - b) / a3;            // real root first
  348.       double real = (((r + s) / -2) - b) / a3;
  349.       double imag = fabs(((r - s) * (sqrt(3.0) / 2)) / a3);
  350.       r2 = CoolComplex(real, imag);        // two complex conjugate
  351.       r3 = CoolComplex(real, -imag);        // roots last.
  352.     } else {                    // resolvent roots are complex
  353.       double rho_3 = curt(rs1.modulus());
  354.       double theta_3 = rs1.argument() / 3.0;
  355.       double rd = 2 * rho_3;
  356.       double cd = ::cos(theta_3) / -2.0;
  357.       double sd = ::sin(theta_3) * sqrt(3.0) / 2.0;
  358.       double a3 = 3 * a;
  359.       if (b < 0) {                // root with largest magnitude
  360.     r1 = CoolComplex(((-2*rd*cd) - b) / a3, 0); // first
  361.     r2 = CoolComplex((rd * (cd + sd) - b) / a3, 0);
  362.     r3 = CoolComplex((rd * (cd - sd) - b) / a3, 0);
  363.       } else {
  364.     r1 = CoolComplex((rd * (cd - sd) - b) / a3, 0);
  365.     r2 = CoolComplex((rd * (cd + sd) - b) / a3, 0);
  366.     r3 = CoolComplex(((-2*rd*cd) - b) / a3, 0);
  367.       }
  368.     }
  369.     return 3;                    // number of roots=3
  370.   }
  371. }
  372.  
  373.  
  374. // roots of quartic -- a*x^4 + b*x^3 + c*x^2 + d*x + e = 0.
  375. //         Return four roots, largest magnitude first.
  376. //         Decompose quartic into two quadratics for better numerical accuracy
  377. //         than directly solving the 4 roots using Ferrari's formula.
  378.  
  379. int CoolComplex::roots_of_quartic (const double& a, const double& b, 
  380.                    const double& c, const double& d, 
  381.                    const double& e,
  382.                    CoolComplex& r1, CoolComplex& r2, 
  383.                    CoolComplex& r3, CoolComplex& r4) {
  384.   if (a < 0) {
  385.     return roots_of_quartic(-a, -b, -c, -d, -e, r1, r2, r3, r4);
  386.   } else if (a == 0) {                // special cases
  387.     return roots_of_cubic(b, c, d, e, r1, r2, r3);
  388.   } else if (e == 0) {
  389.     r4 = 0;
  390.     return roots_of_cubic(a, b, c, d, r1, r2, r3) + 1;
  391.   } else {                    // normal case
  392.     double s;                    // most pos real root of resolvent
  393.     {
  394.       CoolComplex rs1, rs2, rs3;        // roots of resolvent
  395.       roots_of_cubic(1,
  396.              -c,
  397.              (b * d) - (4 * a * e),
  398.              (4 * a * c * e) - ((a * d * d) + (b * b * e)),
  399.              rs1, rs2, rs3);
  400.       if (rs3.imaginary() != 0)            // 2 resolvent roots are imaginary
  401.     s = rs1;
  402.       else if (rs1.real() > rs3.real())        // most positive/negative real 
  403.     s = rs1;                // root is either rs1 or rs3.
  404.       else
  405.     s = rs3;
  406.     }
  407.     double s1 = sqrt((b * b) - (4 * a * (c - s)));
  408.     double s2 = sqrt((s * s) - (4 * a * e));
  409.     double s1s2 = (b * s) - (2 * a * d);
  410.     double a2 = 2 * a;
  411.     if ((s1 * s2 * s1s2) < 0) {            // same sign?
  412.       roots_of_quadratic(a2,
  413.              (b - s1),
  414.              (s + s2),
  415.              r1, r2);
  416.       roots_of_quadratic(a2,
  417.              (b + s1),
  418.              (s - s2),
  419.              r3, r4);
  420.     } else {
  421.       roots_of_quadratic(a2,
  422.              (b - s1),
  423.              (s - s2),
  424.              r1, r2);
  425.       roots_of_quadratic(a2,
  426.              (b + s1),
  427.              (s + s2),
  428.              r3, r4);
  429.     }
  430.     return 4;                    // number of roots=4
  431.   }
  432. }
  433.  
  434.  
  435.